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Abstract 

'We present a massively parallel quantum Monte Carlo based implementation of real-space dynamical mean-field theory for gen- 
eral inhomogeneous correlated fermionic lattice systems. As a first application, we study magnetic order in a binary mixture of 
repulsively interacting fermionic atoms harmonically trapped in an optical lattice. We explore temperature effects and establish 
signatures of the Neel transition in observables directly accessible in cold-atom experiments; entropy estimates are also provided. 
We demonstrate that the local density approximation (LDA) fails for ordered phases. In contrast, a "slab" approximation allows us 
to reach experimental system sizes with Oi lQr') atoms without significant loss of accuracy. 



C/2 



Keywords: dynamical mean-field theory, quantum Monte Carlo, ultracold fermions, optical lattices, antiferromagnetism 
^PACS: 67.85. -d, 03.75. Ss, 71.10.Fd, 71.30.+h 



I 

O 

o 



> 

O 

o 



X 



, Starting with the achievement of Bose Einstein condensation, 
ultracold atomic gases have led to fascinating insights in quan- 
tum many-body phenomena 1 1 ] . With the recent experimental 
successes in realizing quantum degenerate ^ and strongly in- 
teracting Fermi gases |13l|4j] in optical lattices, such systems are 
considered as highly tunable quantum simulators of condensed 
jnatter ijstl • The recent observation of the fermionic Mott tran- 
sition in binary mixtures of "^''/T on cubic lattices Hi marks 
important progress in this respect. The next experimental chal- 
lenge, taken up by many of the leading groups, is to reach and 
identify an ordered antiferromagnetic (AF) Neel phase. 

Ultracold quantum gases on optical lattices have several ad- 
vantages in comparison to solid state systems; however, the de- 
tection methods established for solids do not always find corre- 
spondence in cold-atom experiments. In particular, the detec- 
tion of the AF order parameter is not straightforward. There- 
fore, it is important to identify fingerprints of AF phases that 

directly accessible experimentally. 
■ Due to the intrinsic inhomogeneity of trapped atomic clouds, 
even the qualitative interpretation of experimental data may de- 
pend on corresponding quantitative simulations. The dynamical 
mean-field theory (DMFT) is well established as a powerful, 
nonperturbative approach to interacting Fermi systems 
Spatial inhomogeneities of the optical lattice can be captured 
within a real-space extension of the method (RDMFT) ISHSD or, 
approximately, by applying DMFT within a local density ap- 
proximation (LDA). In either case, the numerical accuracy de- 
pends on the method chosen as DMFT impurity solver. Quan- 
tum Monte Carlo (QMC) based solvers are precise and even 
numerically cheap at or above the temperature ij) ranges rele- 
vant for AF ordering and accessible in cold-atom experiments. 

In this paper, we introduce our Hirsch-Fye QMC IH based 
RDMFT implementation, with focus on computational aspects. 



Then, we turn to a first application, establishing the essential 
physics related to AF order in a relatively small system involv- 
ing (9(100) atoms. In order to minimize finite-size effects, we 
choose a square lattice geometry for which finite-temperature 
antiferromagnetism is not really physical; however, due to the 
mean-field character of the DMFT, the results are representative 
of three-dimensional cubic systems, up to a change in energy 
scales. Using selective simulations with (9(1000) atoms, we 
establish that the properties of large three-dimensional clouds 
can be accurately extracted from simulations for central two- 
dimensional slabs and that simulation boxes can be chosen 
smaller than one would naively expect. 

Model and Methods - Binary mixtures of equivalent 
fermions in an optical lattice are well described by the Hubbard 
model [with (isotropic) trapping potential V,- = Vorf/a^], 
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Here, n,o- = c. c. , c. (c: ) are annihilation (creation) operators 
for a fermion with (pseudo) spin cr e (t, i) at site / (with coordi- 
nates r,), t is the hopping amplitude between nearest-neighbor 
sites (ij), f/ > is the on-site interaction, and fi is the chemical 
potential. We choose Vo - 0.25? and use f as the energy unit. 

The RDMFT approach IH m is based on the following ex- 
pression for the Green function matrix G for Matsubara fre- 
quency ia)„ = (2« + l)i7TT (with spin indices suppressed): 



[G(/w„)],.y = tij + [iaj„ H- yu - y,- - £,('«„)] 6ij ; 



(2) 



here the only approximation is the DMFT assumption of a local 
(i.e. site-diagonal) self-energy S, which corresponds to a mo- 
mentum independence of S in translation-invariant systems (0]. 
The standard DMFT impurity problem, one for each inequiv- 
alent lattice site, is solved in this work using the Hirsch-Fye 
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Figure 1: RDMFT-QMC results for square lattice (Volt = 0.25, Ujt = 10). 
First row: as evidenced by the local magnetization, a large AF core is strongly 
polarized (up to 90%) for T < 0.l4t (left column); with increasing T, both ex- 
tent and magnitude of the AF order decay until a paramagnetic phase is reached 
for r > » 0.32?. Second row: the double occupancy is strongly enhanced 
in AF regions; in contrast, the density (third row) shows little T dependence. 



QMC algorithm OlOll . The discretization of the imaginary-time 
interval < t < 1 /(ksT) introduces a systematic error which 



can be eliminated by extrapolations At ^ 111 U 11211 . In this 
work, we choose Art -0.1 unless indicated otherwise. 

Note that (|2]l only couples G and E at the same Matsubara fre- 
quency; thus, the matrix inversions can be performed in parallel 
for the 500 Matsubara frequencies explicitly taken into account. 
In contrast, the impurity equations couple different frequencies 
(and imaginary times), but are local in the site indices. There- 
fore, this part of the self-consistency problem can be performed 
in parallel for all inequivalent impurities. In addition, the im- 
portance sampling for each impurity may be distributed (via 
MPI) over some 10 - 100 processes. In total, the program can 
saturate 500 CPU cores for large enough problems. 

Results for square lattice - For the square lattice, the nonin- 
teracting dispersion extends from -4f to 4f, leading to a band- 
width of W = 8f. The chosen interaction U = lOt is already 
in the strongly correlated regime: at half filling, the fermions 
would always be localized. However, in the inhomogeneous 
trapped system, an outer shell has to remain delocalized due to 
the outward decay of the filling to zero. Consequently, stag- 
gered magnetic order can only be expected in the core region. 

This is exactly what is seen in Fig.[T](first row): at low T (left 
column), the core shows a nearly perfect staggered magnetiza- 
tion. With increasing T (from left to right), both the extent 
of the ordered region and its polarization decrease until the or- 
der is lost at the Neel temperature Tn ~ 0.32f. Unfortunately, 
this most obvious signature of AF order is not directly acces- 
sible experimentally, primarily due to lack of single-site reso- 
lution. A measurement of particle density profiles (third row) 
also wouldn't help in this case as they are practically unchanged 
(at this scale) through the transition. In contrast, the double oc- 
cupancy, the probability of two particles occupying the same 
lattice site, provides a distinct signal: at high T, it is feature- 
less in the center, with a maximum value of about 0.02. Only 
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Figure 2: Main panel: RDMFT-QMC estimates of entropy for square lattice 
(Vo/t = 0.25, U/t = 10); QMC data for the average entropy per particle (circles, 
for discretization At = 0.1) and for the entropy of the central site (squares, 
also for At = 0.1) are hardly distinguishable from corresponding estimates for 
At = 0.05 (crosses). Inset: temperature derivative dnccnun\/dT of the central 
density for representative values of the chemical potential fi. 



at low temperatures, it is enhanced, by up to 50%, in the de- 
veloping central antiferromagnetic core. For measurements of 
the double occupancy, accuracies of about 0.01 have already 
been established |3]; thus, only a minor refinement is needed to 
resolve this AF signal (when low enough T can be reached). 

Experimentally, the fermionic systems are prepared without 
an optical lattice and with negligible interactions; the initial 
temperature can, therefore, be determined by fitting the density 
profile with a Fermi distribution function. Then, the lattice is 
switched on; due to the localization, this also induces a contact 
interaction. Clearly, the temperature of a final equilibrium state 
will differ greatly from the initial value (in an unknown way); 
however, the process is essentially adiabatic. Therefore, only 
the entropy density is a reliable temperature scale. 

Monte Carlo importance sampling cannot directly access en- 
tropy (or free energy) information. Thus, we revert to the ther- 
modynamic relation dS /dfi = dN/dT and obtain the entropy 
density from integrating over temperature derivatives of parti- 
cle densities (using the fact that fi — » -oo corresponds to an 
empty system with zero entropy). The result is shown in Fig. 
|2] after a steep initial increase, the central entropy density (di- 
amonds) develops a kink around Tn ~ 0.32f (arrow) and then 
reaches a plateau with a value slightly above the mean field 
value of log(2) for a spin system. The average entropy density 
(circles) is always larger due to the low-density shell; in partic- 
ular, the kink is much weaker and a significant slope remains 
above T;^. Corresponding entropy data for smaller discretiza- 
tion (crosses) show excellent agreement, which establishes that 
systematic errors are insignificant. We should point out that fine 
grids in T and fi are essential for this study since strong peaks 
in the integrand dn/dT (see inset of Fig.|2]l have to be resolved. 

One might ask whether expensive RDMFT calculations as 
described above are really necessary for properly resolving the 
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Figure 3: RDMFT-QMC estimates of the AF order parameter, i.e., the stag- 
gered magnetization, for square lattice (U/t = 8, T/t = 0.14, half filling at 
central site) for various trap stengths V (symbols) as a function of the efl'ective 
chemical potential = f^o - Vgrf/a^. Even for nearly realistically weak trap 
strength (Vq = 0.06), the transition is much smoother, by proximity effects, than 
estimated from LDA (using single-site DMFT and QMC; dash-dotted line). 



ordering phenomena. After all, in the case of the paramagnetic 
Mott metal-insulator transition, the LDA, which approximates 
the properties of each site by those of a homogeneous sys- 
tem with the same local potential, yields static observables in 
nearly perfect agreement with RDMFT [9']. However, as seen in 
Fig. |3] the order parameter profiles show enormous deviations 
from LDA even for the largest systems, which fully justifies the 
larger numerical effort of the RDMFT approach. 

Efficient simulations for three dimensions - Let us, finally, 
leave the "toy" case of small two-dimensional systems and dis- 
cuss strategies for realistic simulations of current experimental 
setups involving cubic lattices and 0{\Qr') particles. While such 
system sizes do not pose fundamental difficulties on the impu- 
rity part of the self-consistency problem, due to linear scaling 
and perfect parallelism, a full matrix inversion (with cubic scal- 
ing in the number of sites) as required in (|2]i is out of question, 
at least for dense-matrix algorithms. In our current implementa- 
tion and for typical numerical accuracies, the cost of the matrix 
inversion dominates beyond some 1000 lattice sites. 

Fortunately, the system sizes that have to be simulated ex- 
plicitly can be reduced quite drastically. First of all, periodic 
boundary conditions should be used even for the inhomoge- 
neous trapped systems. Then, the simulation box can be cho- 
sen much smaller than the atomic cloud provided that the AF 
core il3ll does not touch the boundaries, as shown in Fig. |4l 
while spikes are visible (at radial distances of about la) in mjtag 
for a 14 X 14 X 14 system (lower panel, short-dashed lines) at 
Tit - 0.33, the 16 x 16 x 16 data (long-dashed fines) are al- 
ready smooth at this scale, indicating that finite-size effects be- 
come negligible. At T jt - 0.41, closer to the Neel temperature 
Tn ~ 0.45f, the accuracy in mjtag improves further due to the 
smaller core. 

Even larger savings are possible by restricting the simula- 
tion to a central slab, with periodic boundary conditions in the 
perpendicular direction. In practice, we neglect the (small) per- 
pendicular component of the trapping potential in the slab, i.e., 
replace the spherical potential by a cylindrical one, which is 
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Figure 4: RDMFT-QMC estimates for cubic lattice {Volt = 0.05, U/t = 12) 
using periodic boundary conditions, either for a full 3D lattice or using the slab 
approximation (see text). For the density (top panel), the slab results (dashed 
lines) are hardly distinguishable from the full calculations (solid and dotted 
lines); finite-size eft'ects are very small. The agreement is also very good for the 
order parameter (bottom panel), even close to the Neel temperature Tn ~ 0.45?, 
with the exception of spikes corresponding to edges of the 3D systems. Insets: 
cenfi'al density and AF ordering pattern from 22 X 22 slab calculation. 



an excellent approximation for large enough systems. Then, a 
thermodynamic limit is reached, within numerical accuracy, al- 
ready at layer thickness of 4. Corresponding results, indicated 
by solid and dotted lines in Fig.|4] show negligible size depen- 
dence and agree well with the direct 3D results, independent of 
the temperature (density profiles for T/t = 0.41, nearly indistin- 
guishable from those for T/t - 0.33, are omitted). In contrast, 
the LDA (dash-dotted lines) is completely off for mstaa. 

In conclusion, we have implemented the RDMFT approach 
for inhomogeneous correlated Fermi systems on a QMC basis. 
In this work, we have shown, for a relatively small system, that 
the onset of antiferromagnetic order at low T is signaled by an 
enhanced double occupancy. This effect and the impact of the 
interaction strenght U are discussed on a more quantitative level 
for large cubic systems elsewhere [14], using the slab approxi- 
mation established in this paper. 
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